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Abstract 

Motivated by the theory of self-duahty which provides a variational formulation and resolution for non 
self- adjoint partial differential equations [6l|7], we propose new templates for solving large non-symmetric 
linear systems. The method consists of combining a new scheme that simultaneously preconditions and 
symmetrizes the problem, with various well known iterative methods for solving linear and symmetric 
problems. The approach seems to be efficient when dealing with certain ill-conditioned, and highly 
non-symmetric systems. 

1 Introduction and main results 

Many problems in scientific computing lead to systems of linear equations of the form. 

Ax = b where A e M"^" is a nonsingular but sparse matrix, and 6 is a given vector in M", (1) 

and various iterative methods have been developed for a fast and efficient resolution of such systems. The 
Conjugate Gradient Method (CG) which is the oldest and best known of the nonstationary iterative meth- 
ods, is highly effective in solving symmetric positive definite systems. For indefinite matrices, the mini- 
mization feature of CG is no longer an option, but the Minimum Residual (MINRES) and the Symmetric 
LQ (SYMMLQ) methods are often computational alternatives for CG, since they are applicable to systems 
whose coefficient matrices are symmetric but possibly indefinite. 

The case of non-symmetric linear systems is more challenging, and again methods such as CGNE, CGNR, 
GMRES, BiCG, QMR, CGS, and Bi-CGSTAB have been developed to deal with these situations (see the 
survey books [9] and [llj). One approach to deal with the non-symmetric case, consists of reducing the 
problem to a symmetric one to which one can apply the above mentioned schemes. The one that is normally 
used consists of simply applying CG to the normal equations 

A^Ax = A^b or AA^y = b, x = A^y. (2) 

It is easy to understand and code this approach, and the CGNE and CGNR methods are based on this idea. 
However, the convergence analysis of these methods depends closely on the condition number of the matrix 
under study. For a general matrix A, the condition number is defined as 

n{A)^\\A\\-\\A-% (3) 
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and in the case where A is positive definite and symmetric, the condition number is then equal to 



HA) = (4) 

where Amin(^) (resp., Amax(^)) is the smallest (resp., largest) eigenvalue of A). The two expressions can 
be very different for non-symmetric matrices, and these are precisely the systems that seem to be the most 
pathological from the numerical point of view. Going back to the crudely symmetrized system ([2]), we echo 
Greenbaum's statement [9] that numerical analysts cringe at the thought of solving these normal equations 
because the condition number (see below) of the new matrix A is the square of the condition number of 
the original matrix A. 

In this paper, we shall follow a similar approach that consists of symmetrizing the problem so as to be 
able to apply CG, MINRES, or SYMMLQ. However, we argue that for a large class of non-symmetric, 
ill-conditionned matrices, it is sometimes beneficial to replace problem ([T]) by one of the form 

A^MAx = A^Mb, (5) 

where M is a symmetric and positive definite matrix that can be chosen properly so as to obtain good 
convergence behavior for CG when it is applied to the resulting symmetric A'^ MA. This reformulation 
should not only be seen as a symmetrization, but also as preconditioning procedure. While it is difficult to 
obtain general conditions on M that ensure higher efficiency by minimizing the condition number k{A'^ MA), 
we shall show theoretically and numerically that by choosing M to be either the inverse of the symmetric 
part of A, or its resolvent, one can get surprisingly good numerical schemes to solve (H]). 
The basis of our approach originates from the selfdual variational principle developed in [SI [7] to provide a 
variational formulation and resolution for non self-adjoint partial differential equations that do not normally 
fit in the standard Euler-Lagrangian theory. Applied to the Unear system ([T]) , the new principle yields the 
following procedure. Split the matrix A into its symmetric Aa (resp., anti-symmetric part Aa) 

A = As + Aa, (6) 

where 

As ■^1{A + A^) and A^ := ^(^ - ^^)- (7) 

Proposition 1.1 (Selfdual symmetrization) Assume the matrix A is positive definite, i.e., for some S > 0, 

{Ax,x) >S\x\^ for all X £R". (8) 

The convex continuous functional 

I{x) = ^{Ax,x) + ^{A-\b- AaX),b- AaX) - {b,x) (9) 
then attains its minimum at some x in R", in such a way that 

I{x) = inf I{x) = (10) 
Ax = b. (11) 

Symmetrization and preconditioning via selfduality: Note that the functional / can be written as 

I{x) = \{Ax,x) + {AaA-H^b,x) + \{A-\b), (12) 

where 

A:^As- AaA-^Aa = A^A~^A. (13) 
By writing that DI{x) — 0, one gets the following equivalent way of solving ([1]). 
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// both A G R"^" and its symmetric part As are nonsingular, then x is a solution of the equation if and 
only if it is a solution of the linear symmetric equation 

A^A-^Ax = {As - AaA;^Aa)x ^ b - AaA^^b = A^A^^b. (14) 

One can therefore apply to ((M]) all known iterative methods for symmetric systems to solve the non- 
symmetric linear system ([T]). As mentioned before, the new equation (jl4p can be seen as a new symmetriza- 
tion of problem ([1]) which also preserves positivity, i.e., A^A~^A is positive definite if A is. This will then 
allow for the use of the Conjugate Gradient Method (CG) for the functional /. More important and less 
obvious than the symmetrization effect of A, is our observation that for a large class of matrices, the con- 
vergence analysis on the system (ITi|) is often more favorable than the original one. The Conjugate Gradient 
method -which can now be applied to the symmetrized matrix A— has the potential of providing an efficient 
algorithm for resolving non-symmetric linear systems. We shall call this scheme the Self-Dual Conjugate 
Gradient for Non-symmetric matrices and we will refer to it as SD-CGN. 

As mentioned above, the convergence analysis of this method depends closely on the condition number k{A) 
oi A = A^Aj^A which in this case is equal to k{A). We observe in section 2.3 that even though k{A) could be 
as large as the square of k{As), it is still much smaller that the condition number of the original matrix k{A). 
In other words, the inverse C of A^A^^ can be an efficient preconditioning matrix, in spite of the additional 
cost involved in finding the inverse of As. Moreover, the efficiency of C seems to surprisingly improve in 
many cases as the norm of the anti-symmetric part gets larger (Proposition 2.2). A typical example is when 
the anti-symmetric matrix Aa is a multiple of the symplectic matrix J (i.e. J J* = — = /). Consider then 
a matrix A^ = As -\- which has an arbitrarily large anti-symmetric part. One can show that 

^i{A,)<K{As)-\-e^XrnUAf, (15) 

which means that the larger the anti-symmetric part, the more efficient is our proposed selfdual precondi- 
tioning. Needless to say that this method is of practical interest only when the equation AsX = d can be 
solved with less computational effort than the original system, which is not always the case. 
Now the relevance of this approach stems from the fact that conjugate gradient methods for nonsymmetric 
systems are costly since they require the storage of previously calculated vectors. It is however worth noting 
that Concus and Golub [3 and Widlund ,15j have also proposed another way to combine CG with a pre- 
conditioning using the symmetric part As-, which does not need this extended storage. Their method has 
essentially the same cost per iteration as the preconditioning with the inverse of A^ Aj^ that we propose for 
SD-CGN and both schemes converge to the solution in at most N iterations. 

Iterated preconditioning: Another way to see the relevance of As as a preconditioner, is by noting that 
the convergence of "simple iteration" 

AsXk = -AaXk-i +b (16) 

applied to the decomposition of A into its symmetric and anti-symmetric parts, requires that the spectral 
radius p{I — Aj-^A) — p{Aj^Aa) < 1. By multiplying by A~^ , we see that this is equivalent to the 
process of applying simple iteration to the original system ([T|) conditioned by Aj^, i.e., to the system 

A:' Ax = A;'b. (17) 

On the other hand, "simple iteration" applied to the decomposition of A into As and AaA^^Aa is given by 

AsXk = AaA-^AaXk-i + b - AaA-^b. (18) 

Its convergence is controlled by p{I — A~^A) ~ p{{Aj^ Aa)"^) = p{Aj-^ Aa)"^ which is strictly less than 
p{Aj^ Aa) , i.e., an improvement when the latter is strictly less than one, which the mode in which we have 
convergence. In other words, the linear system (jl4[) can still be preconditioned one more time as follows: 

// both A e R"^" and its symmetric part As are nonsingular, then x is a solution of the equation |ip if and 
only if it is a solution of the linear symmetric equation 

Ax := A-^A^A-^Ax = [/ - {As^Aaf]x = {I - A;^Aa)A;H = A^^A'^A^^b. (19) 
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Note however that with this last formulation, one has to deal with the potential loss of positivity for the 
matrix A. 

Anti-symmetry in transport problems: Numerical experiments on standard linear ODEs (Example 3.1) 
and PDEs (Example 3.2), show the efhciency of SD-CGN for non-selfadjoint equations. Roughly specking, 
discretization of differential equations normally leads to a symmetric component coming from the Laplace 
operator, while the discretization of the non-self-adjoint part leads to the anti-symmetric part of the coeffi- 
cient matrix. As such, the symmetric part of the matrix is of order O(^), while the anti-symmetric part is 
of order O(^), where h is the step size. The coefficient matrix A in the original system ([1]) is therefore an 
0{h) perturbation of its symmetric part. However, for the new system (1141) we have roughly 

A^As- AaA-'Aa = O(-l) - 0{\)0{e)0{^) = 0{^) - 0(1), (20) 

making the matrix A an 0(1) perturbation of As, and therefore a matrix of the form Ag + al becomes a 
natural candidate to precondition the new system (|14p . 

Resolvents of Ag as preconditioners: One may therefore consider preconditioned equations of the form 
A^MAx = A'^Mb, where M is of the form 

Af„ - {aA, + (1 - a)l) or Np = l3Aj^ + (1 - /?)/, (21) 

for some < a, /? € M, and where / is the unit matrix. 

Note that we obviously recover ^ when a = 0, and HH) when a = 1. As a — > the matrix aAg + (1 — a)! 
becomes easier to invert, but the matrix 

Ai^c, = A^{aAs + {l-a)I)''^A (22) 

may become more ill conditioned, eventually leading (for a = 0) to A^ Ax = A^h. There is therefore a 
trade-off between the efficiency of CG for the system ([5]) and the condition number of the inner matrix 
aAs -|- (1 — a)/, and so by an appropriate choice of the parameter a we may minimize the cost of finding a 
solution for the system ([T]). In the case where As is positive definite, one can choose -and it is sometimes 
preferable as shown in example (3.4)- a > 1, as long as a < jup — , where A^j^^ is the smallest eigenvalue 
of As- Moreover, in the case where the matrix A is not positive definite or if its symmetric part is not 
invertible, one may take a small enough, so that the matrix Ma (and hence ^i,q) becomes positive definite, 
and therefore making CG applicable (See example 3.4). Similarly, the matrix Np = I3A~'^ + (1 ~ P)I provides 
another choice for the matrix M in ([5]), for /3 < ^^"""^-^ where A^g^^ is the largest eigenvalue of As- Again 
we may choose a close to zero to make the matrix Njs positive definite. As we will see in the last section, 
appropriate choices of /?, can lead to better convergence of CG for equation (O . 
One can also combine both effects by considering matrices of the form 

La,i3 = {aAs + (1 - a)l) + [31, (23) 

as is done in example (3.4). 

We also note that the matrices M'^ := {aA's + (1 - a)iy^ and N'p := /?(A^)"^ + (1 - /?)/ can be other 
options for the matrix Af , where A'^ is a suitable approximation of As, chosen is such a way that M'^q and 
N'^q can be relatively easier to compute for any given vector q. 

Finally, we observe that the above reasoning applies to any decomposition A — B + C oi the non-singular 
matrix A £ M"^", where B and (B — C) are both invertible. In this case, B{B — C)^^ can be a preconditioner 
for the equation ([T]). Indeed, since B — CB~^C = {B ~ C)B^^A, a; is a solution of ([T]) if and only of it is a 
solution of the system 

[B - C)B-^Ax = {B - CB-^C)x = b- CB-^b. (24) 

In the next section, we shall describe a general framework based on the ideas explained above for the use 
of iterative methods for solving non-symmetric linear systems. In section 3 we present various numerical 
experiments to test the effectiveness of the proposed methods. 
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2 Selfdual methods for non-symmetric systems 



By selfdual methods we mean the ones that consist of first associating to problem ([T]) the equivalent system 
([5]) with appropriate choices of M, then exploiting the symmetry of the new system by using the various 
existing iterative methods for symmetric systems such as CG, MINRES, and SYMMLQ, leading eventually 
to the solution of the original problem In the case where the matrix M is positive definite, one can then 
use CG on the equivalent system ([S]). This scheme (SD-CGN) is illustrated in Table (1) below, in the case 
where the matrix M is chosen to be the inverse of the symmetric part of A. If M is not positive definite, then 
one can use MINRES (or SYMMLQ) to solve the system (HH). We will then refer to them as SD-MINRESN 
(i.e., Self-Dual MINRES for Nonsymmetric linear equations). 



2.1 Exact methods 

In each iteration of CG, MINRES, or SYMMLQ, one needs to compute Mq for certain vectors q. Since 
selfdual methods call for a conditioning matrix M that involves inverting another one, the computation of 
Mq can therefore be costly, and therefore not necessarily efficient for all linear equations. But as we will 
see in section 3, M can sometimes be chosen so that computing Mq is much easier than solving the original 
equation itself. This is the case for example when the symmetric part is either diagonal or tri-diagonal, 
or when we are dealing with several linear systems all having the same symmetric part, but with different 
anti-symmetric components. Moreover, one need not find the whole matrix M, in order to compute Mq. 
The following scheme illustrates the exact SD-CGN method applied in the case where the coefficient matrix 
A in IT]) is positive definite, and when {As)~^ Aq can be computed exactly for any given vector q. 



Given an initial guess xq, 
Solve AsU = b 
Compute b ^ b — AaU- 
Solve AsUo = AaXQ 

Compute tq — b — AsXq + AaUo and set po — tq. 
For k=l,2, . . . , 
Solve AsZ = AaPk-i 
Compute w = AsPk^i — AaZ . 

Set Xk = Xk-i + ak-iPk-i, where Uk^i = <p^_i.t„> 
Cpmpute Tfe = rk-i - ak-iw. 

Set pf,^rk+ bk-iPk-i, where bk-i = <rfcliirt^i> ■ 
Check convergence; continue if necessary. 



<rfc_i ,rfc_i> 



Table 1: GCGN 

In the case where A is not positive definite, or when it is preferable to choose a non-positive definite 
conditioning matrix M, then one can apply MINRES or SYMMLQ to the equivalent system ([5]). These 
schemes will be then called SD-MINRESN and SD-SYMMLQN respectively. 

2.2 Inexact Methods 

The SD-CGN, SD-MINRESN and SD-SYMMLQN are of practical interest when for example, the equation 

AsX = q (25) 

can be solved with less computational effort than the original equation ([T]). Actually, one can use CG, 
MINRES, or SYMMLQ to solve ^ in every iteration of SD-CGN, SD-MINRESN, or SD-SYMMLQN. 
But since each sub- iteration may lead to an error in the computation of (|25p . one needs to control such 
errors, in order for the method to lead to a solution of the system H]) with the desired tolerance. This 



5 



leads to the Inexact SD-CGN, SD-MINRESN and SD-SYMMLQN methods (denoted below by ISD-CGN, 
ISD-MINRESN and ISD-SYMMLQN respectively). 

The following proposition -which is a direct consequence of Theorem 4.4.3 in 9:- shows that if we solve 
the inner equations ([25]) "accurately enough" then ISD-CGN and ISD-MINRESN can be used to solve (P) 
with a pre-determined accuracy. Indeed, given e > 0, we assume that in each iteration of ISD-CGN or 
ISD-MINRESN, we can solve the inner equation -corresponding to Ag- accurately enough in such a way 
that 

UA,--AaA;^Aa)p-iA,p~Aay)\\ - \\AaA;^AaP-Aay\\ <e, (26) 
where y is the (inexact) solution of the equation 

AsV = Aap. (27) 

In other words, we assume CG and MINRES are implemented on in a finite precision arithmetic with 
machine precision e. Set 

.o:=2(. + 4). .,:=2(7 + nMl_W^)^^ (28) 

where \D\ denotes the matrix whose terms are the absolute values of the corresponding terms in the matrix D. 
Let Ai < ... < A„ be the eigenvalues of {As — AaAj-^Aa) and let Tk+i^k be the (fc -I- 1) x fc tridiagonal matrix 
generated by a finite precision Lanczos computation. Suppose that there exists a symmetric tridiagonal 
matrix T, with Tk+i^k as its upper left {k + 1) x k block, whose eigenvalues all lie in the intervals 

^ = ur=i[A,-5,A, + (5], (29) 

where none of the intervals contain the origin, let d denote the distance from the origin to the set S, and let 
Pk denote a polynomial of degree k. 



Proposition 2.1 The ISD-MINRESN residual r{^-^ then satisfies 

(30) 



< V(l + 2eo)(A: + l) min max ^^(0)! + 2Vfc(-f )ei. 

ro Pk z=s d 



If A is positive definite, then the ISD-CGN residual r^'^ satisfies 

^-jr-TT ^ V(l + 2eo)(A„+,5)/d min max |pfc(z)| + Vfc(-f )ei. (31) 
I Foil d 

It is shown by Greenbaum [B] that Tk+i^k can be extended to a larger symmetric tridiagonal matrix T whose 
eigenvalues all lie in tiny intervals about the eigenvalues of {Ag — AaAj^Aa). Hence the above proposition 
guarantees that if we solve the inner equations accurate enough, then ISD-CGN and ISD-MINRESN con- 
verges to the solution of the system [T] with the desired relative residual (see the last section for numerical 
experiments) . 



2.3 Preconditioning 

As mentioned in the introduction, the convergence of iterative methods depends heavily on the spectral 
properties of the coefficient matrix. Preconditioning techniques attempt to transform the linear system ([1]) 
into an equivalent one of the form C^^Ax — C^^b, in such a way that it has the same solution, but hopefully 
with more favorable spectral properties. As such the reformulation of (1) as 

A^Aj^Ax = A^A-^b, (32) 

can be seen as a preconditioning procedure with C being the inverse of A'^A^^. The spectral radius, and 
more importantly the condition number of the coefficient matrix in linear systems, are crucial parameters for 
the convergence of iterative methods. The following simple proposition gives upper bounds on the condition 
number of ^ = A^Aj^A. 
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Proposition 2.2 Assume A is an invertible positive definite matrix, then 

< min{Ki, K2}, (33) 



where 



Proof: We have 



We also have 



:= k(^) + JfJ^ and := + ^^^II^aI) ■ (34) 



^ ^niax{As^ 



Since n{A) — ^"""=(4) ^ it follows that k{A) < ni. 
To obtain the second estimate, observe that 



Amin(^) — ^min{As — A^A^ Aa) > Amm(^s) 



- x^Ax X*{As — AaA ^Aa)x 

XmaA^) - sup— ^= sup r—^ 

x^O \X\ x^Q Fl 



^min {Ag ) 



i(^) — ^mi7i{As ~ AaAg ^ Aa) > ^min{~AaAg ^ Aa) 

■ r '^^'^ -^aAg ^ AaX 
— mt = 

x^O X X 

.^^^, {AaX)'^A^\AaX) ^ {AaX)^{AaX) 

™o'^ {AaX)^{AaX) ^ X'^X 
^ . ^ {AaX)^ Aj\AaX) . ^ xT'{Aa)'^{Aa)x 

> mt — -— — TTfrr-. — ^ — X mi = 

x^O (AaX)^ [AaX) x=iO X^ X 

X '\min{{Aa)'^ Aa) 



^jnax {A 



\ ( A ) 

With the same estimate for Xmaxi-A-) we get k{A) < K2 



Remark 2.1 Inequality ((33)) shows that SD-CGN and SD-MINRES can be very efficient schemes for a large 
class of ill conditioned non-symmetric matrices, even those that are almost singular and with arbitrary large 
condition numbers. It suffices that either ki or K2 be small. Indeed, 

• The inequality k{A) < ki shows that the condition number k{A) is reasonable as long as the anti- 
symmetric part Aa is not too large. On the other hand, even if \\Aa\\ is of the order of Ai„ax(^s), and 
k{A) is then as large as k{As)'^, it may still be an improved situation, since this can happen for cases 
when k{A) is exceedingly large. This can be seen in example 2.2 below. 

• The inequality k{A) < K2 is even more interesting especially in situations when X^in{— A"^) is arbitrarily 
large while remaining of the same order as This means that k{A) can remain of the same order 
as k{As) regardless how large is Aa- 

A typical example is when the anti-symmetric matrix Aa is a multiple of the symplectic matrix J (i.e. 
J J* = —J^ = I). Consider then a matrix A^ — As + which has an arbitrarily large anti-symmetric 
part. By using that k{A) < K2, one gets 

k{A,) < k{As) + e2A„,ax(As)'. (35) 

Here are other examples where the larger the condition number of A is, the more efficient is the proposed 
selfdual preconditioning. 
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Example 2.2 Consider the matrix 



A, = 



-1 

-l + e 



(36) 



which is a typical example of an ill-conditioned non-symmetric matrix. One can actually show that k{A^) ~ 
0(i) ^ cxD as e with respect to any norm. However, the condition number of the associated selfdual 
coefficient matrix 

1 \ 

Ae = As AaiAs) Aa = 

[ e 

is k{A^) = and therefore goes to 1 as e ^ 0. Note also that the condition number of the symmetric 
part of A^ goes to one as e — > 0. In other words, the more ill-conditioned problem IH) is, the more efficient 
the selfdual conditioned system dH]) is. 

We also observe that k{Aj^A) goes to oo as e goes to zero, which means that besides making the problem 
symmetric, our proposed conditioned matrix Aj^A has a much smaller condition number than the matrix 
A^^A, which uses As as a preconditioner. 

Similarly, consider the non-symmetric linear system with coefficient matrix 



A, = 



-I + e 
-1 



(37) 



As e — !■ 0, the matrix becomes again more and more ill-conditioned, while the condition number of its 
symmetric part converges to one. Observe now that the condition number of A^ also converges to 1 as e 
goes to zero. This example shows that self-doual preconditioning can also be very efficient for non-positive 
definite problems. 



3 Numerical Experiments 

In this section we present some numerical examples to illustrate the proposed schemes and to compare them 
to other known iterative methods for non-symmetric linear systems. Our experiments have been carried out 
on Matlab (7.0.1.24704 (R14) Service Pack 1). In all cases the iteration was started with xq = 0. 

Example 3.1 Consider the ordinary differential equation 

~ey" + y' = f{x), on [0,1], 2/(0) = 2/(1) = 0. (38) 

By discretizing this equation with stepsize 1/65 and by using backward difference for the first order term, 
one obtains a nonsymmetric system of linear equations with 64 unknowns. We present in Table 2 below, the 
number of iterations needed for various decreasing values of the residual e. We use ESD-CGN and ISD-CGN 
(with relative residual 10~^ for the solutions of the inner equations) . We then compare them to the known 
methods CGNE, BiCG, QMR, GGS, and BiGGSTAB for solving non-symmetric linear systems. We also 
test preconditioned version of these methods by using the symmetric part of the corresponding matrix as a 
preconditioner. 
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Table 2: Number of iterations to find a solution with relative residual 10 ^ for equation ([38)1 . f{x) is chosen 
so that y — X sin(7ra;) is a solution. 



N=64 


e = 10-" 


e = 10-^ 


e = 10-* 


e = 10-" 


e = 10-^*^ 


€ = 10-^" 


ESD-CGN 


22 


8 


5 


4 


3 


2 


ISD-CGNflO-') 


24 


9 


6 


4 


3 


2 


GCNE 


88 


64 


64 


64 


64 


64 


QMR 


114 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


PQMR 


34 


51 


50 


52 


52 


52 


BiCGSTAB 


63.5 


78.5 


92.5 


98.5 


100.5 


103.5 


PBiCGSTAB 


26.5 


46.5 


50.5 


50 


51.5 


51.5 


BiCG 


125 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


PBiCG 


31 


44 


50 


50 


52 


52 


CGS 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


PCGS 


27 


51 


46 


46 


46 


48 



Table 3: Number of iterations to find a solution with relative residual 10 ^ for equation ((38)l . f{x) is chosen 
so that y — ^ggT^? is a solution, while the stepsize used is 1/129. 



N=128 


e = 10-" 


e = 10-^ 


e = 10-* 


e = 10"" 


e = lO-i'^ 


e= 10-^" 


ESD-CGN 


37 


11 


6 


4 


3 


2 


ISD-CGN(10-0 


38 


12 


7 


4 


3 


2 


GCNE 


266 


140 


128 


128 


128 


128 


QMR 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


PQMR 


40 


77 


87 


92 


90 


85 


BiCGSTAB 


136.5 


167.5 


241 


226.5 


233.5 


237.5 


PBiCGSTAB 


35.5 


87.5 


106.5 


109 


110.5 


110.5 


BiCG 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


PBiCG 


37 


76 


84 


89 


85 


91 


CGS 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


> 1000 


PCGS 


34 


80 


96 


91 


94 


90 



As we see in Tables 2 and and 3, a phenomenon similar to Example \2.S\ is occuring. As the problem gets 
harder (e smaller), SD-CGN becomes more efficient. These results can be compared with the number of 
iterations that the HSS iteration method needs to solve equation 138\) (Tables 3,4, and 5 in 12]). 

Example 3.2 Consider the partial differential equation 

Oil 

-Au + aix,y)— = f{x,y), < x < I, < y < 1, (39) 
with Dirichlet boundary condition. 

The number of iterations that ESD-CGN and ISD-CGN needed to find a solution with relative residual 10-", 
are presented in Table 4 below for different coefficients a(x, y). 
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Table 4: Number of iterations (I) for the backward scheme method to find a solution with relative residual 
10-6 for equation ^ (Example 3.2) 



a(x,y) 


N 


I (ESD-CGN) 


I (ISD-CGN) 


Solution 


100 


49 


18 


18 


random 


100 


225 


40 


37 


random 


100 


961 


44 


46 


random 


100 


961 


52 


51 


sin7ra;sin7r';/.expffa;/2 + w)'^) 


1000 


49 


10 


10 


random 


1000 


225 


31 


31 


random 


1000 


961 


36 


37 


random 


1000 


961 


31 


39 


sin7rxsin7r2/.exp((a;/2 + y)'^) 


10*^ 


49 


4 


4 


random 


10*^ 


225 


6 


6 


random 


10« 


961 


6 


6 


random 


10« 


961 


6 


6 


sin7r2;sin7r?/.exp((a;/2 + y)"^) 


10^« 


961 


2 


2 


sin7rxsin7r2/.exp((a;/2 + y)"^) 



Table 5: Number of iterations (I) for the centered difference scheme method for equation ([39| (Example 3.2) 



a(x,y) 


N 


I (ESD-CGN) 


Solution 


Relative Residoual 


1 


49 


21 


random 


6.71 X 10"*^ 


1 


225 


73 


random 


9.95 X 10-*^ 


1 


961 


91 


random 


8.09 X 10-*^ 


1 


961 


72 


sin7r2;sin7r?/.exp((a;/2 + y)'^) 


9.70 X 10-" 


10 


49 


18 


random 


9.97 X 10-" 


10 


225 


65 


random 


5.90 X 10-" 


10 


961 


78 


random 


8.95 X 10-" 


10 


961 


65 


sinTTCCsinTT?/. exp((a;/2 + y)'^) 


7.78 X 10-" 


100 


49 


31 


random 


6.07 X 10-" 


100 


225 


42 


random 


5.20 X 10-" 


100 


961 


43 


random 


5.03 X 10-" 


100 


961 


38 


sin7rxsin7r?/.exp((a;/2 + y)'^) 


4.69 X 10-" 


1000 


49 


65 


random 


4.54 X 10-" 


1000 


225 


130 


random 


8.66 X 10-" 


1000 


961 


140 


random 


2.12 X 10-" 


100 


961 


150 


sin TTX sin iry. exp((a;/2 + y)"*) 


5.98 X 10-" 



Table 4 and 5 can be compared with Table 1 in [15 , where Widlund had tested his Lanczos method for 
non-symmetric linear systems. Comparing Table 5 with Table 1 in [15] we see that for small a{x,y) (1 and 
10) Widlund's method is more efficient than SD-CGN, but for large values of a, SD-CGN turns out to be 
more efficient than Widlund's Lanczos method. 

Remark 3.3 As we see in Tables 2,3, and 4, the number of iterations for ESD-CGN and ISD-CGN (with 
relative residual 10-^ for the solutions of the inner equations) are almost the same One might choose dynamic 
relative residuals for the solutions of inner equations to decrease the average cost per iterations of ISD-CGN. 
It is interesting to figure out whether there is a procedure to determine the accuracy of solutions for the inner 
equations to minimize the total cost of finding a solution. 
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Example 3.4 Consider the partial differential equation 

-Au + lo ^ieM^M^' + y'h) ^ io^^p(3 5(^2 ^ y2y^du ^ [0, 1] x [0, 1], (40) 

ox ox 

with Dirichlet boundary condition, and choose f so that sm^nx) sin^iry) e'xp{{x/2 + y)^) is the solution of the 
equation. We take the stepsize h = 1/31 which leads to a linear system Ax = b with 900 unknowns. Table 5 
includes the number of iterations which CG needs to converge to a solution with relative residual 10^® when 
applied to the preconditioned matrix 

A^{aA-^ + {l~a)I)A. (41) 

Table 5 can be compared with Table 1 in '151, where Widlund has presented the number of iterations needed 
to solve equation |^Q[ ). 



Table 6: Number of iterations for a solution with relative residual 10 ^ for example 3.3 when SD-CGN is 
used with the preconditioner ([1T|) for different values of a. 





I 




\s / l-a \ 


I 


oo{a — 0) 


> 5000 




0.1 


232 


0(a = 1) 


229 




0.2 


237 


-0.1 


221 




0.4 


249 


-0.25 


216 




0.8 


263 


-0.5 


201 




1 


272 


-0.7 


191 




5 


384 


-0.8 


186 




10 


474 


-0.9 


180 




20 


642 


-0.95 


179 




50 


890 


-0.99 


177 




100 


1170 


-0.999 


180 




1000 


2790 


-0.9999 


234 




10000 


4807 



Remark 3.5 As we see in Table 5, for ^max{'^~^) — ~-99 "we have the minimum number of iterations. 
Actually, this is the case in some other experiments, but for many other system the minimum number of 
iterations accrues for some other a with —1 < ^max{^~^) — 0- experiments show that for a well chosen 
a > 1, one may considerably decrease the number of iterations. Obtaining theoretical results on how to choose 
parameter a in \41\ seems to be an interesting problem. 

Note that the coefficient matrix of the linear system corresponding to (HO]) is positive definite. Hence we 
may also apply CG with the preconditioned symmetric system of equations 

A^iA, ~ aA^,i„/)-M = A^iAs - aX^^.jr'b, (42) 

where X^^^ is the smallest eigenvalue of Ag and a < 1. The number of iterations function of a, that CG 
needs to converges to a solution with relative residual 10^^ are presented in Table 7. 
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Table 7: Number of iterations to find a solution with relative residual 10 ^ for equation pO]) when SD-CGN 
is used with the preconditioner (j42|) for different values of a. 



ex. 


I 





229 


0.5 


204 


0.9 


177 


0.99 


166 


0.999 


168 


0.9999 


181 


0.99999 


194 


0.999999 


222 


0.9999999 


248 


0.99999999 


257 



Remark 3.6 As we see in the above table, for a = 0.99 in we have the minimum number of iterations. 
Obtaining theoretical results on how to choose the parameter a seems to be an interesting problem to study. 

We also repeat the experiment by applying CG to the system of equations 

(a, - 0.99X:^,J)-' - ^l) A = A^ ({A, - o.99A^,„/)"i - ^l) b. (43) 

Then CG needs 131 iterations to converge to a solution with relative residual 10~^. 
As another experiment we apply CG to the preconditioned linear system 

A:^A^A;^A = A;^A^A;\ 

to solve the non-symmetric linear system obtained from discritization of the Equation (|40p . The CG converges 
in 31 iterations to a solution with relative residual less than 10^^. Since, we need to solve two equations with 
the coefficient matrix As, the cost of each iteration in this case is towice as much as SD-CGN. So, by the 
above preconditioning we decrease cost of finding a solution to less that 62/131 of that of SD-CGN (System 

Example 3.7 Consider now the following equation 

-Au + io ^(^^P(3-5(^' + ^')") ^ i0exp(3.5(:r^ + y^))^ ~ 200^ ^ /(x), on [0, 1] x [0, 1], (44) 
ox ox 

If we discretize this equation with stepsize 1/31 and use backward differences for the first order term, we 
get a linear system of equations Ax = b with A being a non-symmetric and non-positive definite coefficient 
matrix. We then apply CG to the following preconditioned, symmetrized and positive definite matrix 

A'' {{A, - oAf„i„/)-i + (3I)A - A^((^ - aA^i„/)"i + /3I)b, (45) 

with a < 1. For different values of a the number of iterations which CC needs to converge to a solution with 
the relative residual 10^^ are presented in Table 8. 



12 



Table 8: Number of iterations to find a solution with relative residual 10 ^ for equation (|44]) when SD-CGN 
is used with the preconditioner (j45|) for different values of a and (3. 



a 


^ = 


B -.99/A" 

r 1 max 


10 


543 


424 


5 


446 


352 


2.5 


369 


288 


1.5 


342 


264 


1.1 


331 


258 


1.01 


327 


259 


1.001 


333 


271 


1.0001 


368 


289 


1.00001 


401 


317 



We repeat our experiment with stepsize 1/61 and get a system with 3600 unknowns. With a = —1.00000001 
and (3 — 0, CG converges in one single iteration to a solution with relative residual less than 10~^. We also 
apply QMR, BiCGSTAB, BiCG, and CGS (also preconditioned with the symmetric part as well) to solve the 
corresponding system of linear equations with stepsize 1/31. The number of iterations needed to converge 
to a solution with relative residual 10^^ are presented in Table 9. 

Table 9: Number of iterations to find a solution with relative residual 10^^ for equation (j44|l using various 
algorithms. 



N=900 


I 


CGNE 


> 5000 


QMR 


3544 


PQMR 


490 


BiCGSTAB 


> 5000 


PBiCGSTAB 


Breaks down 


BiCG 


4527 


PBiCG 


> 1000 


CGS 


1915 


PCGS 


649 
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